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I .  INTRODUCTION 


Since  an  ELF  signal  from  a  remote  transmitter  is  received  over  a  range  of 
azimuth  angles/  lateral  ionospheric  gradients  produced  by  sporadic  E  layering 
or  nuclear  depressions  can  produce  significant  effects  on  propagation  in  the 
lower  ELF  band.  This  is  because  the  Fresnel  zone  size  can  be  large  as 
compared  with  the  distance  over  which  the  ionosphere  changes  significantly  in 
the  lateral  direction.  Although  a  number  of  workers  have  addressed  the 

12  3 

question  of  off-path  effects  (Wait  ,  Galejs  ,  Greifinger  and  Greifinger  , 
Field  Field  and  Joiner  8,7 ,  Pappert  8)  no  formulation  exists  which  can 
fully  account  for  the  propagation  effects  produced  by  a  localized  disturbance 
with  simultaneous  allowance  for  vertical  inhomogeneity,  lateral  inhomogeneity, 
and  anisotropy  in  a  spherical  geometry.  It  has  been  common  practice  to 
estimate  the  effects  of  lateral  gradients  by  using  a  simple  surface  propaga¬ 
tion  model  introduced  by  Wait  and  more  fully  developed  by  the  Greifingers  and 
Field.  The  formulation  reduces  the  problem  to  an  integral  equation  descrip¬ 
tion  of  propagation  along  the  earth's  surface.  The  theory  is  predicated  on 
the  palatable  assumption  that  the  field  can  be  separated  into  lateral  and 
height  dependent  functions  when  the  lateral  ionspheric  gradients  are  consider¬ 
ably  smaller  than  the  vertical  gradients.  When  applying  the  method  to 
nocturnal  environments  additional  assumptions  are  made.  Among  these  is  the 
omission  of  nonreciprocal  effects.  This  is  well  justified  in  the  ambient 
case9  as  well  as  for  daytime  and  depressed  ionospheres.  However,  it  is  known 
that  under  sporadic  E  layering  considerable  mixing  between  TE  and  TM  wave  can 
occur Thus,  when  the  surface  propagation  model  is  applied  to  sporadic  E 
environments  the  scattered  TE  component  is  neglected.  The  validity  conditions 
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for  the  formulation  are  probably  best  satisfied  under  conditions  of  either 
natural  or  man-made  depressed  ionospheres. 

The  purpose  of  this  report  is  to  document  a  computer  program,  based  on 
the  surface  propagation  model,  which  is  useful  to  the  user  community  for 
estimating  the  effects  of  localized  ionospheric  disturbances  on  propagation  of 
the  vertical  electric  field  component  Ez  at  lower  ELF  frequencies.  A  moments 
method11  serves  as  the  basis  for  the  solution  of  the  integral  equation. 
Though  the  method  is  powerful,  in  the  present  program,  practical  storage 
requirements  restrict  application  to  disturbances  which  effectively  vanish 
outside  a  rectangle  of  several  megameters  on  one  side.  It  is  hoped  that 
limitation  will  be  relaxed  in  future  work.  The  program  allows  in  some 
measure,  for  modeling  of  rectangular,  circular,  and  elliptical  disturbance 
shapes.  The  lateral  propagation  function,  W,  defined  as  the  ratio  of  the 
disturbed  laterally  dependent  part  of  the  vertical  field  component,  Ez,  to  the 
undisturbed  field  component  is  calculated  as  is  the  absolute  value  of  the 
total  Ez  field  component  in  the  disturbed  guide.  The  latter  calculation 
allows  in  approximate  manner  for  guide  height  effects  via  WKB  formalism  as 
applied  to  waveguide  propagation. 

The  program  requires  eigenangle  inputs  for  both  the  ambient  and  disturbed 
regions  of  the  guide  as  well  as  end-on  horizontal  dipole  excitation  factors 
for  the  vertical  Ez  field  component.  These  must  be  supplied  from  a  waveguide 
program  such  as  that  of  reference  12. 
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II.  SUMMARY  OF  EQUATIONS 


Subject  to  the  assumption  that  the  vertical,  Ez»  field  component  can  be 
separated  into  lateral  and  height  dependent  functions,  the  lateral 
dependence  <p(x,y)  is  given  by^ 

2  00 

<(>(x,y)  *  ^(x.y)  -  /  /  dx'dy'  (s2{x',y')  -  S2)  G(|r-r'  |)  ♦(x*  ,y' )  M) 

«CO 

where  S  is  the  sine  of  the  eigenangle  for  the  disturbed  guide,  SQ  the  sine  for 
the  unperturbed  guide,  k  the  free  space  wave  number,  and  the  superscript  i 
signifies  the  unperturbed  incident  field.  The  Green's  function  G(|r-r' |)  is 
given  by 


where 

r  =  xi  +  yj  and  r*  *  x'i  +  y* j.  (3) 

Here  a@  is  the  earth's  radius  and  the  x,  y,  z  are  the  Cartesian  coordinates,  z 
is  measured  vertically  upwards  into  the  ionosphere  and  x  is  measured 
horizontally  with  x-z  the  plane  of  incidence.  Unit  vectors  in  the  x  and  y 
directions  are  denoted  by  1  and  j  .  The  square  root  factor  in  equation  (2) 
has  been  introduced  to  allow  for  the  geometric  spreading  appropriate  to  a 
spherical  geometry.  Beyond  that,  x,  y  and  x' ,  y'  are  taken  to  be  rectangular 
coordinates  in  a  flat  earth  geometry.  The  quantity  Hp(2)  tfte  Hankel 
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coordinates  in  a  flat  earth  geometry.  The  quantity  2  ^  is  the  Hankel 
function  of  order  p  of  the  second  kind.  The  unperturbed  lateral  dependent 
function,  ,  is  taken  to  be  appropriate  for  the  vertical  electric  field,  Ez, 
launched  by  the  en<2hon  component  of  a  ground  based  horizontal  electric  dipole 
source  oriented  in  the  x  direction  and  given  by 


-5 


a  sin 
e  v  a 


(¥) 


(4) 


Ir| 


where  B  is  a  constant.  Again  the  square  root  factor  allows  for  the  spreading 
appropriate  to  a  spherical  geometry.  The  discrete  analog  of  equation  (1)  for 
field  points  within  the  disturbed  region  consists  of  the  N  x  N  system  of 
linear  equations: 


N 

l  '  ®  =  1*2,..., N  (5) 

,  mn  n  m 
n=l 


Allowance  is  made  for  modeling  rectangular,  circular,  and  elliptical  distur¬ 
bances  by  subdividing  them  into  square  mesh  cells.  It  is  common  practice  in 
such  circumstances  to  simplify  the  integrations  involving  cylindrical 
functions  by  approximating  each  square  cell  by  a  circle  of  equal  area.  It  is 
also  common  practice  to  take  the  electric  field  to  be  constant  over  the  area 
of  a  cell.  However,  larger  cells  can  be  tolerated  if  the  electric  field  is 
allowed  to  vary  over  the  area  of  a  cell.  At  least  two  methods  of  allowing  for 
this  variation  have  been  described  in  the  literature13.  In  this  work  the 
method  called  'plane  wave  correction'  is  used.  The  method  assumes  isotropic 
propagation  within  and  between  each  cell  and  that  should  be  a  reasonable 
approximation  in  the  present  problem  since  anisotropy  of  propagation  at  ELF 
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frequencies  is  quite  small.  Translated  into  the  notation  of  this  work,  the 
results  of  reference  13  yield  for  the  A-matrix  elements: 


-f  4 

The  Green's  function,  G,  is  given  by  equation  (2)  with  r  replacing  r  and  r 

m  n 

replacing  r'  .  Also,  Jp  is  the  Bessel  function  of  order  p  of  the  first 

kind.  In  terms  of  the  *¥  determined  by  equation  (5),  and  the  A^  given  by 

equation  (7),  the  field  at  a  point,  r  ,  exterior  to  the  disturbed  region  is 

ra 

given  by: 


N 

-  1 

n=1 


A 

mn  n 


(8) 


Equations  (5)  through  (8)  are  used  in  the  present  program  to  determine 
the  dB  value. 


W  -  20  1°910I'Hx,y)/^i(x,y)  }, 


<9) 
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of  the  disturbed  lateral  function  relative  to  its  undisturbed  value. 

Another  output  of  the  program  makes  allowance  for  dependence  of  the 
vertical  electric  field  on  height  of  the  guide  via  the  approximate  WKB 
formalism.  For  a  laterally  homogeneous  guide  the  Ez  field  generated  by  the 
end-on  component  of  a  horizontal  dipole  may  be  expressed  as 


E 

z 


_  s3/2  (1-jRjjRJ 

~  ”lQ  3f  z 

URl 


{c(i- Vix.y) 


(10) 


where  F  is  the  modal  function  and  9F/90  its  derivative  evaluated  at  the  eigen- 
angle.  S  and  C  are  the  sine  and  cosine  of  the  eigenangle.  ^R^  and  n  R ^  are 
TE  and  TM  Fresnel  reflection  coefficients  referenced  to  the  ground  and  ^R^  is 
the  ionospheric  TE  reflection  coefficient  referenced  to  the  ground.  In  the 
absence  of  anisotropy  the  quantity  (1-jR  ^R^)  would  cancel  an  identical  tei.n 
occurring  in  the  (9F/90)  .  The  Fresnel  coefficient  ^Ru  is 


where  is  the  complex  refractive  index  of  the  ground.  Because  the  magnitude 
of  Nq  is  much  greater  than  unity  in  the  lower  ELF  band  good  approximations 
are: 


1  +  jR,  -  2  and  0(1-^,)  =  2/Ng 


(12) 


Thus,  within  the  spirit  of  the  WKB  approximation  the  Ez  field  becomes 
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Ez  -  ~4Qi 


[  3/2  _  1  1/2  f  3/2  _  1/2 

3f  ^  ~1R11RJ  I  I  3f  ^  "lRilR.L^  IT 

36  Jr  L  36  JtGt 


»(x,y)  (13) 


where  Q  is  a  constant  dependent  upon  dipole  moment  and  frequency  and  the 

subscripts  r  and  t  stand  for  receiver  and  transmitter.  That  is,  the  first 

term  in  parenthesis  is  evaluated  at  the  receiver  while  the  second  term  in 

parenthesis  and  the  factor  N  ^  are  evaluated  at  the  transmitter. 

G 

The  factor  B  in  equation  (4)  is  taken  to  be 


*kS0ae  -(3/4) ui 
- -  e 


To  express  the  vertical  electric  field,  E^ ,  as  given  by  equation  (13)  in 
microvolts/m,  the  factor  Q  has  the  value 


Q  =  2.849  x  10'3  f3''2  (Idl) 
kHz 


where  f^^  is  the  frequency  in  kHz  and  Idl  is  the  current  moment  in  ampere 


meters. 
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III.  GEOMETRICAL  MODELING  OF  THE  DISTURBANCE 


As  described  in  the  introduction,  the  program  can  be  used  to  model 
square,  rectangular,  circular,  and  elliptical  disturbances.  In  all  cases  the 
disturbance  is  defined  to  be  symmetrical  about  both  its  x  and  y  axes.  The 
special  case  of  a  rectangular  disturbance  is  required  to  be  uniformly  dis¬ 
turbed.  Let  Xq  and  yp  denote  the  coordinates  of  the  center  of  the  disturbance 
with  respect  to  the  location  of  the  transmitter.  Let  1^  and  Ly  be  the  size  of 
the  disturbance  along  the  x  and  y  axis.  The  disturbance  is  overlaid  by  a  grid 
which  has  n}{  squares  along  the  x  axis  and  ny  squares  along  the  y  axis.  The 
choice  of  Lx,  Ly,  nx,  and  ny  must  be  such  that  Lx/nx  =  Ly/ny. 

Since  we  assume  that  the  disturbance  is  symmetrical  about  the  x  and  y 
axis,  we  only  need  to  specify  the  waveguide  eigen  solution  parameters  along 
the  x  axis.  Let  us  denote  these  solutions  by  Si.  We  take  i  =  0  to  denote  the 
ambient  or  undisturbed  values.  The  remaining  (S.,  through  S^)  is  assumed 
to  be  uniformly  distributed  along  the  x  axis  from  x^  to  Xq  +  V2  Lx«  Note  that 
the  value  of  nx  is  not  related  to  that  of  N.  If  N  =  1,  then  a  uniform  distur¬ 
bance  will  be  assumed.  This  is  required  for  the  rectangular  disturbance.  The 
remaining  problem  is  to  fill  the  disturbance  grid  with  interpolated  values  of 
S. 

Let  us  now  consider  a  single  subsquare  within  a  square  disturbance.  Let 
the  coordinates  of  the  center  of  this  subsquare  be  x  and  y.  The  smaller  value 
of  |x  -  xQ|  and  ly  -  yQ|  is  used  to  interpolate  a  value  of  S  from  the  input 
list  of  S^.  This  results  in  the  disturbance  grid  being  filled  as  illustrated 
in  figure  1.  In  the  figure  the  similarly  shaded  regions  would  all  have  the 
same  value  of  S. 
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Figure  1 .  Diagram  illustrating  the  distribution  of  S  within  a  rectangular  disturbance. 
The  similarly  shaded  regions  have  the  same  value  of  S. 


Let  us  consider  an  elliptical  disturbance  with  an  axial  ratio  R  = 

L  /L  .  The  center  of  a  given  subsquare  is  at  x  and  y.  We  can  define  an 

y  x 

ellipse  which  is  concentric  with  that  of  the  outer  edge  of  the  disturbance. 
This  ellipse  intersects  the  x  axis  at  a  distance  from  Xg  defined  by 


{(x-x  )2+  [ <y~yn)/R]2  K2 


The  above  expression  can  be  used  to  calculate  the  value  of  A  for  each  corner 
of  the  subsquare,  say  A^  where  i  =  1,  2,  3  or  4.  If  the  smallest  of  the  A^  is 
greater  than  Lx,  then  the  subsquare  is  entirely  outside  the  disturbance  and  SQ 
is  used  in  the  grid  at  that  subsquare.  If  the  largest  of  the  A^  is  less  than 
1<X,  then  the  subsquare  is  entirely  inside  the  disturbance  and  the  value  of  S 
for  that  subsquare  is  interpolated  from  the  input  list  of  S^.  The  remaining 
case  is  that  of  a  subsquare  on  the  edge  of  the  disturbance.  In  this  case  the 
subsquare  is  subdivided  into  16  smaller  squares.  Let  the  coordinates  of  the 
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For  each  of  these 


center  of  each  of  these  smaller  squares  be  x_  and  y_ . 

m  m 

squares  we  calculate  ^  using  equation  (16).  The  value  of  S  for  the  subsquare 
is  given  by 

nS  -  ( 16-n)Sft 


where  n  is  the  number  of  subsquares  which  are  within  the  ellipse.  A  circular 
disturbance  is  treated  the  some  way  as  an  elliptical  one  with  R  =  1 . 
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IV.  DESCRIPTION  OF  INPUT 


All  input  to  the  sporadic-E  program  is  given  in  a  data  deck  on  the 
standard  input  unit.  A  listing  of  sample  input  showing  the  data  deck  setup  is 
shown  in  figure  2. 

There  are  two  parts  to  the  input.  This  first  part  is  read  in  by  means  of 
a  FORTRAN  NAMELIST  input  format.  The  first  card  of  each  set  of  input  must 
contain  '&DATUM'  in  columns  2-7.  This  is  followed  by  at  least  one  blank  and 
then  data  items  separated  by  commas.  The  data  items  have  the  following  forms: 
'variable  name'  =  constant, 
or 

'array  name'  =  set  of  constants  (all  separated  by  commas). 

The  data  list  is  terminated  by  'SEND'.  All  cards  must  have  a  blank  in  column  1. 

The  second  part  of  the  input  follows  the  NAMELIST.  The  first  card  for 
this  part  is  an  identification  card.  It  contains  up  to  80  columns  of  alpha¬ 
numeric  information  and  is  used  to  label  the  output  plots.  Following  the 
identification  card  a  series  of  punched  cards  (obtained  from  the  programs 
described  in  reference  14  with  NPUNCH  =  1 )  is  input  for  each  mode. 

The  first  card  gives  the  values  of  R,  FREQ,  AZIM,  CODIP,  MAGFLD,  SIGMA 
and  EPSR.  Next  there  are  two  cards  per  mode.  The  first  of  these  contains  the 
complex  eigenangle  at  the  ground  and  values  for  the  complex  quantities  T1  and 
T2.  The  second  card  contains  the  eigenangle  at  the  ground  (duplicate  input) 
and  values  for  T3  and  T4.  The  quantities  T1,  T2,  T3,  and  T4  are  defined  in 
reference  14. 
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The  following  variables  and  arrays  may  be  specified  in  the  NAMELIST 


input: 

DM 

OMIN 

DMAX 

DELD 

YMAX 

DELY 

I  FLAG 

IGRID 

XO  ,  YO 

NUMX , 
NUMY 

SIZEX 

SIZEY 

I  PLOT 


dipole  moment  in  ampere  meters. 

the  minimum  range  in  kilometers  at  which  fields  are 
calculated,  printed,  and  plotted  for  I FLAG  =  2. 

the  maximum  range  in  kilometers  at  which  fields  are  calculated, 
printed,  and  plotted  for  I FLAG  =  2. 

the  increment,  in  kilometers  at  which  fields  are  calculated, 
printed,  and  plotted  for  IFLAG  =  2. 

the  maximum  off  axis  value,  in  kilometers,  at  which  fields  are 
calculated,  printed,  and  plotted  for  IFLAG  =  1 . 

the  increment  in  kilometers,  at  which  fields  are  calculated 
printed  and  plotted  for  IFLAG  =  1 . 

for  IFLAG  =  1  the  program  calculates  fields  as  a  function  of  y 
-  the  distance  is  fixed  at  DMIN  and  the  disturbance  moves  in  the  y 
direction. 

for  IFLAG  =  2  the  program  calculates  fields  as  a  function  of 
distance  with  the  disturbance  fixed. 

IGRID  =  0  indicates  that  the  disturbance  is  either  square 
or  rectangular. 

IGRID  =  1  indicates  that  the  disturbance  is  either  circular  or 
elliptical. 

coordinates  in  kilometers  at  center  of  disturbance.  YO  is 

also  the  initial  value  at  which  fields  are  calculated,  printed, 

and  plotted  for  IFLAG  =  1 . 

the  disturbance  is  divided  into  NUMX  grids  in  the  x-  direction 
and  NUMY  grids  in  the  y-direction  (these  are  n  ,  n  in  the  text). 

x  y 

defines  the  physical  size  of  the  disturbance.  It  is  SIZEX 
kilometers  by  SIZEY  kilometers. 

a  flag  controlling  whether  or  not  plots  are  generated.  If 
IPLOT  »  0  no  plots  are  generated.  If  IPLOT  =  1  two  plots  are 
generated.  For  IFLAG*2  the  first  plot  is  WMAG(DB)  vs  X(KM) , 
equation  (9).  The  second  plot  consists  of  two  curves:  EZUMAG 
(DB),  equations  (13)  and  (4),  is  a  solid  curve  and  EZPMAG  (DB), 
equations  (13)  and  (8),  is  a  dashed  curve.  For  IFLAG-1  WMAG, 
EZUMAG,  and  EZPMAG  are  plotted  with  respect  to  y  for  a  fixed  value 
of  x. 
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XLNG  the  length  in  inches  £or  the  x-axis  and  y-axis  for  the 
YLNG  field  plots. 

HMIN  the  minimum  and  maximum  values,  in  dB,  desired  on  y-axis  for 
UMAX  the  WMAG  plot. 

EMIN  the  minimum  and  maximum  values  in  DB  desired  on  y-axis  for 

EMAX  the  EZREL  plot. 

XYIC  the  units  per  tic  mark  along  the  x-axis  and  y-axis  for 
YTIC  both  plots. 


Initial  values  of  the  namelist  variables  are  presented  in  table  1. 


TABLE  1 

Namelist  variables  and  initial  values. 


NAME 

INITIAL  VALUE 

UNITS 

DM 

6.75E6 

Ampere  meters 

DMIN 

25.0 

Kilometers 

DMAX 

1000.0 

Kilometers 

DELD 

25.0 

Kilometers 

YMAX 

5000.0 

Kilometers 

DELY 

25.0 

Kilometers 

I  FLAG 

2 

IGRID 

0 

XO 

0.0 

Kilometers 

YO 

0.0 

Kilometers 

NUM 

0 

NUMY 

0 

SIZEX 

1000.0 

Kilometers 

SIZEY 

1000.0 

Kilometers 

I  PLOT 

0 

XLNG 

5.0 

Inches 

YLNG 

6.0 

Inches 

WMIN 

-6.0 

Decibels 

UMAX 

2.0 

Decibels 

EMIN 

-60 

Decibels 

EMAX 

2.0 

Decibels 

XTIC 

200.0 

Kilometers 

YTIC 

0.2 

Decibels 

16 


V.  DESCRIPTION  OF  OUTPUT 


A  listing  of  sample  output  is  shown  in  figure  3.  The  resulting  plots  are 
found  in  figures  4  through  7.  Figures  4  and  5  are  the  output  from  IFLAG  =  2. 

Figures  6  and  7  are  the  output  from  IFLAG  =  1.  The  first  section  of  output  is 

an  echoing  of  the  namelist  input  variables.  This  is  followed  by  a  schematic 
drawing  of  the  disturbed  region  showing  grid  numbers  and  the  coordinates  of 
the  center  of  the  corner  meshes.  The  next  section  shows  additional  input 
parameters:  the  identification  label  that  will  be  on  the  plots,  the 

frequency,  conductivity,  dielectric  function,  and  the  complex  eigenangle  and 
excitation  factor  for  each  region.  YMID  represents  the  y-value  coordinate  of 
the  midpoint  of  the  disturbance. 

In  the  printout,  the  tables  for  WI,  EZO,  and  EZS,  are  the  magnitudes  of 

the  quantities  as  given  by  equations  9,  10,  and  13.  These  quantities  are 

computed  at  the  midpoints  of  each  grid  square. 

The  last  table  in  the  printout  lists  the  following  quantities  at  the 
given  distance  from  the  transmitter  along  the  x-axis  at  y*0: 

EZREL , EZANG  magnitude ( dB)  and  phase  angle  (radians)  of  equation  8 

W MAG, WANG  magnitude( dB)  and  phase  angle  (radians)  of  equation  9 

EZUMAG , EZUANG  magnitude ( dB)  and  phase  angle  (radians)  of  equation  10 

EZPMAG , EZPANG  magnitude( dB)  and  phase  angle  (radians)  of  equation  13 
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Figure  3.  Continued. 
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Figure  3.  Continued. 
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WMAG  (dB) 
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dB  ABOVE  1  MV/m 


SAMPLE  X-VARIATION  PLOT 
SOUFIRE  OR  RECTANGULAR  DISTURBANCE 
FREQ  -  0.075 

SIZEX-  500.0  XO-  2000.0  NUMX-  4 
SIZEY-  1000.0  YO-O.O  NUMY-  8 

Figure  5,  Sample  x-variation  plot  fot  EZUMAG  and  EZPMAG  assuming  a  rectangular  disturbance. 
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WMAG  (dB) 
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SAMPLE  Y-VARIATION  PLOT 
SQUARE  OR  RECTANGULAR  DISTURBANCE 
FREQ  -  0.075  DMIN  -  2000.0 
SIZEX-  1000.0  XO-  1000.0  NUMX-  5 
SIZEY-  1000.0  YO-  0.0  NUMY-  5 


Figure  6.  Sample  y-variation  plot  for  WMAG  assuming  a  rectangular  disturbance. 
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SAMPLE  Y-VARIATION  PLOT 
SQUARE  OR  RECTANGULAR  DISTURBANCE 
FREQ  -  0.075  DMIN  -  2000.0 
SIZEX-  1000.0  X0-  1000.0  NUMX-  5 
SIZEY-  1000.0  Y0-  0.0  NUMY-  5 


Figure  7.  Sample  y-variation  plot  for  EZUMAG  and  EZPMAG  assuming  a  rectangular  disturbance. 
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VI .  PROGRAM  CHECKS 


Several  program  checks  have  been  made  of  flat  earth  geometry  cases  which 
can  be  solved  in  terms  of  well  known  functions.  In  each  case  the  disturbance 
is  azimuthally  symmetric  with  the  transmitter  located  at  the  origin. 

The  first  case  considered  is  that  of  a  uniform  circular  disturbance  for 

which 


_2  2 

S  =  S 
P 

2  2 

s  =so 


for  r<r„ 


for  r>r„ 


(18) 


For  this  case  the  solution  for  the  ratio  of  the  disturbed  lateral  function, 
¥  ,  to  the  undisturbed,  ¥^  is  for  r<r 


h!2V( 


7T 


<kSpr)  - 


(19) 


lHi2)(kspV  h(ksor0)  -  «Vo>  Hj2)(ksoro 
(Ji|kspro>  h(ks0r0)  “  3(kspro)  Hi2>(ksoro^ 


j^(kSr) 


and  for  r> r 


„f  “i2>(>Vol  i^Vo1  -  hlkVo>  I 

°1  »!2’<«oco>  J«“pro>  -  wkVo>  J1,kVo(  I 


(20) 
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In  these  equations 


h(x)  [H<2)(x,  -  H*2)(x)]  (21) 

j(x)  =  j  [Jq (x)  -  J2(x) ]  (22) 

It  is  clear  from  equation  (20)  that  is  constant  for  r^r^  as  expected. 
Figure  8  shows  the  results  calculated  by  using  equation  (19)  along  with  the 
moment  method  results.  The  radius  r^  is  500  km.  The  unperturbed  eigenangle 
is  (83 .985° ,-34 .909° )  or  equivalently  Sq  *  (1 . 185 ,-0 .0681i)  and  the  perturbed 
eigenangle  is  (59.393° ,-65.552® )  or  equivalently  Sp  = (1. 488,-0. 718i)  .  Sq  is 
appropriate  to  a  nighttime  ambient  ionosphere  at  75  Hz,  whereas,  Sp  is  appro- 

Q 

priate  to  a  nighttime  ionosphere  with  a  sporadic  E  layer  .  The  13  x  13  mesh 
results  give  agreement  to  within  a  few  hundreths  of  a  dB  of  the  exact 
results.  In  the  moments  method,  the  lateral  function  is  calculated  at  the 
center  points  of  each  square  mesh  and  the  lateral  function  is  linearly 
interpolated  between  those  points.  Because  of  approximations  made  in  the  slab 
containing  the  transmitter,  the  first  meaningful  data  point  obtained  from  the 
moment  method  is  the  first  point  to  fall  in  a  slab  adjacent  to  that  containing 
the  transmitter.  This  explains  the  starting  ranges  for  the  moment  method 
results. 

A  second  check  case  considered  is  that  of  a  circular  disturbance  which  is 
uniform  out  to  a  radius  rQ,  then  falls  off  as  1/r2  between  rQ,  and  r^»  and  is 
equal  to  SQ  beyond  r^ .  The  mathematical  description  of  S  is 
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UNIFORM  CIRCULAR  DISTURBANCE  (500km  RADIUS  -  XMTR  AT  CENTER! 


Figure  8.  Comparison  between  analytic  solution  and  the  computer  program  output  for  problem  1 :  a  uniform  circular  disturbance 


2  2 
S  “  S 
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2  2 

S  +  (S 

0  p 

■J 


,  r  <  r0 

S0  >  r0  /r'  J  r0  <  r  <  rl 
;  r  >  r. 


When  the  transmitter  is  at  the  center  of  the  disturbance  the  solution  for  the 
ratio  of  the  disturbed  lateral  function,  ^  ,  to  the  undisturbed,  'll1  ,  is  for 


A.  _ 

i 

S0Hj2l(kS0r) 

1 

soHS2)(ksor) 

JL, 

T__ 

so 


In  these  equations 


fsH*2)(kS  r)  +  aJ„(kS  r) ]  ;  r  < 
L  1  p  1  P 


[bJ  (kS.rl  +  dJ  (kS  rl]  ;  r  <  r  <  r 
L  0  J  ~vv  O'1  0  1 


;  r  <  r. 


V=  [1  -  (s2  -  S2)  (kr0)2]V2 


SH<12  ^Spr0]  Jv^kS0r0^  J-v^S0r0^ 

Sh1<kVo)  VkSoV  ^-v^kS0r0) 

0  Jv^kSorl)  J-v^kS0r^ 

0  MkSOrl)  i-v^kS0r1^ 


a  “  A 


-»;2)(ksori) 

“hl(kS0r1) 


I 


Kj|X 


Vk)  =  !  ^Hq-1 (x)  “  (34) 

Again  it  is  clear  from  equation  (26)  that  ty/il'*  is  constant  for  r>r1  .  Figure 
9  shows  results  for  rQ  =  300  km,  r1  =  1000  km  and  the  same  SQ  and  Sp  used  for 
the  results  of  figure  8.  The  moment  method  results  are  for  a  10  x  10  mesh  and 
the  results  beyond  100  km  agree  with  the  exact  values  to  better  than  a  tenth 
of  a  dB. 

A  third  check  based  on  the  model  described  by  equations  (23)  is  shown  in 
figure  10.  SQ  and  are  assigned  the  same  values  as  for  figure  8;  however, 
Tq  has  been  taken  equal  to  588  km  and  r^  =  2000  km.  The  moment  method  result 
is  for  a  17  x  17  mesh  which  is  the  largest  the  program  can  handle  because  of 
storage  limitations.  The  mesh  size  in  this  instance  is  approximately  235  km. 

In  this  connection  Hagmann  et  al13<  state  that  the  mesh  size  must  be  less  than 
155  Ag  where  A^  is  the  unperturbed  wavelength.  In  the  present  case  of  75  Hz, 
Ag=  4000/SQr  =  3375  km  .  This  would  give  a  maximum  cell  size  of  =  523  km. 
It  has  been  our  experience  that  approaching  this  limit  leads  to  substantial 
error  and  it  would  probably  be  best  at  75  Hz  not  to  exceed  mesh  sizes  of 
several  hundred  kilometers.  This  would  limit  the  linear  dimension  of  the 
disturbance  to  something  on  the  order  of  5000  km  at  75  Hz.  Figure  10  shows 
the  agreement  between  the  exact  calculation  and  the  moment  method  to  be  within 
a  few  tenths  of  a  dB.  It  is  also  very  likely  true  that  the  cases  checked  are 
some  of  the  most  difficult  for  the  moment  method  to  handle  because  the 
gradient  of  the  incident  field  is  largest  close  to  the  transmitter.  Thus,  it 
would  be  expected  that  for  disturbances  remote  from  the  transmitter  better 
accuracy  would  be  obtained  for  the  same  mesh  size. 
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VIII.  APPENDIX:  PRCX3RAM  LISTING 
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bHECLQma  ptas  bunk-hot  nim> 


*>FTN,USFO 
FTN  1  OR  1  A 


SPORAO .MAIN 
09/22/82-1  1 : 32 ( 9 , 1 0  ) 

1  .  C 

2.  C  SPCPAPIC-E  program  for  elf 

3  .  C 

4.  IMPLICIT  COMPLEX  (A-H.O-Zt 

5.  C 


6  . 

7  .  C 

8  .  C 
9 .  C 

10 .  C 

11.  C 
12  . 

13  . 

14  . 

15  . 

16  .  C 

17  . 

18  . 

19  . 

20  . 

21  . 

22  . 

23  . 

24  . 

25  . 

26  . 

27  . 

28  .  C 


PARAMETER  MAX NX=12, NDIM -100, NRPTS=500 

MAXNX  IS  The  MAX  NUMBER  OF  DISTURBED  EIGENS  WHICH  CAN  BE  INPUT 

NDIM  IS  THE  MAX  VALUE  OF  NUMX*NU'-’Y 

NRPTS  IS  THE  MAX  NUV3ER  OF  POINTS  TO  PLOT 


COMPLEX  IM/  I  0  .  EO  ,  1  .  EO  )/  ,  J1  .  J2  ,  NGMM  ,  KSO  ,NGSO, 

$  AINDIM.NDIM)  . 

$  £20 ( NDIM) , EZS(NOIM) ,SD<  VDIW) , XTRD<  NOIM) , WI (NQIM) , 

$  SI( MAXNX) .XTRI(MAXNX) 


REAL 

$ 

$ 

S 

$ 

$ 

$ 

s 

$ 

s 

$ 


XM. IV, XGSO .AI.SIZEX.SIZEY.XDINC.XO.YO.VOVERX.DELTAX, 
DElTAY.X1.Y1.XINC.VINC.EPSR, RHOMN, EZREL . EZANG . X A , DE  L  Y , 
YMAX . F , WM AG. WANG . YMID . OMEGA .RADIUS.WAVENO, FREO .SIGMA . 

E  RP-  .  DIST  .DM  IN,  DM.  DELD.  DMAX  .  H  X  I  NC  .  H  Y  1 NC  ,  RG  .  A  IMA  X  ,  DST  , 
RGMAX .RGMIN.EZUMAG. EZUANG. EZPMAG.EZPANG, PI , DTR . EPSO, 

XG(  4  I  . YG( 4  |  . 

X(NDIM) , Y( NDIM) ,CLNO( ND IM ) . X X ( NDIM ) , Y Y ( ND IM ) , 

XI ( MAXNX) , Y I (MAXNX) , 

XP LOT (NRPTS ) ,Y1 PLOT (NRPTS) . Y 2 P LOT ( NRPTS )  , Y3PL0T ( NRPT  S ) , 

ext ic.eytic.wxtic.wytic.xmax.xlng, ylng, wmin.wmax.emin, 

EV.  v 


29.  INTEGER  IROW(NDIM) 

30  .  C 


31.  CHARACTER’1  LABELO/'?'/ 

32.  character*-;  i o ; 20 ) 

33.  CHAP  C  TER*  5  L A  BE  L 1 / '  ?'/.  LABEL2/' _ ? ' / . LABE L3/ '  __  '/ 

34.  CHAR AC  TE  R  *  8  XLABEL 

35.  CHAR  AC  T  E  R  *  26  LABEL 

36.  C 

37  .  C0I.1-.10N/0NE/H12  ,H22,J1  .  J2  .  SO  .  ARGO  ,  K  A  .  KSO 

38  .  COMMON,  I-'  CELS/  ID.  LA3EL.X  LABEL 


39  .  CO'.  '.'UN/  1 N  KU  i  ,  FREC.DMIN,  DMAX  .SIZE-'  .  SI  ZE  Y  ,  XO  ,  YO  ,  F.XTIC  ,  EYTIC  ,  WXTIC  . 

40.  $  WYTIC.XLNG.YLNG.EMIN.EMAa.WMIN.WMAX.XMAX.NUMX.NUMY, 

41  .  S  I  FLAG 


42 .  COMMON/ PL DA f A/ X P LOT , Y 1  PLOT , Y2P LOT , Y3P LO T , NPTS 

43  .  C 

44.  C  DM  -  OIPGLE  MOMENT  IN  AMPERE-METERS 

45.  C  DM! N  -  MINIMUM  RANGE  (KM)  FOR  DISTANCE  VARIATION 

46.  C  DMAX  -  MAXIMUM  RANGE  (KM)  FOR  DISTANCE  VARIATION 

47.  C  DELD  -  INCREMENT  ( KV I  FOR  DISTANCE  VARIATION 

48.  C  Yf/.AX  -  MAXIU.,M  OFF-AXIS  VALUE  (KM)  FOR  Y  VARIATION 

49.  C  DELY  -  INCREMENT  (KM)  FOR  Y  VARIATION 

50.  C  I F  L  AG= 1  Y  VARIATION 

51.  C  I F  L AG  =  2  DIST  VARIATION 

52.  C  IGRID  =  0  SOUARE  OR  RECTANGULAR  DISTURBANCE 

53.  C  IGRID  =  1  CIRCULAR  OR  ELLIPTICAL  DISTURBANCE 

54.  C  XO.YO  -  COORDINATES  (KM)  AT  CENTER  OF  DISTURBANCE 

55.  C  YO  IS  ALSO  INITIAL  VALUE  FOR  Y  VARIATION 
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1 

1 

1 


t 

1 


56.  C 

57.  C 

58.  C 

59.  C 

60 .  C 
61  .  C 

62 .  C 

63.  C 
64  .  C 

65.  C 

66.  C 

67 .  C 

68 .  C 

69.  C 

70  . 

71  . 

72  . 

73.  C 

74. 

75  . 

76. 

77  . 

78. 

79. 

80  . 

81  . 

82. 

83. 

84 .  C 

85.  C 

86. 

87.  C 
88  . 

89. 

90. 

91  .  C 

92. 

93. 

94. 

95. 

96. 

97  . 

98  . 

99. 

100  . 

101  .  C 
102 .  C 

103  . 

104  . 

105. 

106. 

107 .  C 

108.  C 

109. 

110. 

Ill  . 

112. 


NUMX, NUMY  -  THE  DISTURBANCE  IS  DIVIDED  INTO  NUMX  GRIDS  IN  THE 
X  DIRECTION  AND  NUMY  GRIDS  IN  THE  Y  DIRECTION 
SI2EX , SI2EY  -  DEFINES  THE  PHYSICAL  SI2E  OF  THE  DISTURBANCE  -  IT  IS 
SI2EX  KM  BY  S12EY  KM 
IPlOT=0  DON'T  PLOT 
IPLOT  =  1  PLOT 

XLNG.YLNG  -  LENGTH  (INCHES)  OF  X-AXIS  AND  Y-AXIS  RESPECTIVELY 
WMIN.WMAX  -  MINIMUM  AND  MAXIMUM  (DB)  FOR  Y-AXIS  OF  WMAG  PLOT 
EMIN.EMAX  -  MINIMUM  AND  MAXIMUM  (DB)  FOR  Y-AXIS  OF  E2  PLOT 
EXTIC.EYTIC  -  UNITS  PER  TIC  MARK  ALONG  X  AND  Y  AXIS  RESPECTIVELY 
FOR  E2  PLOTS 

WXTIC.WYTIC  -  UNITS  PER  TIC  MARK  ALONG  X  AND  Y  AXIS  RESPECTIVELY 
FOR  WMAG  PLOTS 

NAME  LI  ST/ DATUM/DM, DMI N. DMAX ,DELD, YMAX, DELY , IF  LAG, I  GRID, 

$  XO, YO, NUMX, NUMY, S 2 2EX.S12EY, IPLOT, XLNG.YLNG, 

$  WMIN,WMAX,EMIN,EMAX.EXTIC,EYTIC,WXTIC,WYTIC 


data  iplot/o/ 

DATA  PI/3 . 14159265358 979E0/.DTR/. 0174532 92 EO/ , EPSO/8 .85434E-1 2/ 
DATA  XLNG/5. /, YLNG/6. / . WMI N/-6 . / . WMAX/2 . / , EMI N/-6 . / , EMAX/2 . / 
DATA  EXT1C/200./.EYTIC/5./.SI2EX/1 .E3/.SI2EY/1 .E3/ 

DATA  WXT1C/200./.WYT1C/1 ./ 

DATA  DMIN, 25 . E 0/ . DMAX/ 1 000 . EO/ , DE LD/25 . EO/ 

DATA  YMAX/SOO. EO/ , DELY/25 . EO/ . XO/0 . EO/ , YO/O.EO/ 

DATA  DM/6.75E6/ 

DATA  IFLAG/2/, IGRID/O/ 

DATA  CONS T/(-. 707106781 I86548E0,-. 7071 06781 186548E0)/ 


DEFINE  CAPG2(ARG)  »  P1/(2.0E0*  SIN( ARG/6371 ,0E0) ) 

PRINT  915 
READ(5, DATUM) 

WR I T  E ( 6 , DA  TUM ) 

NU  *  NUMX  *NUMY 
I F ( NU  ,GT.  NDIM) 

S  THEN 

PRINT  971 
STOP 

END  IF 
NPTS  *  0 
DIST  *  DMIN 
YMID  =  YO 

SET  UP  DISTURBED  REGION  GRID 

XINC  *  S I 2EX/NUMX 

HXINC  *  . 5*X INC 

YINC  *  SI2EY/NUMY 

HYINC  *  . 5* Y I NC 

GRID  ELEMENTS  MUST  BE  SQUARE 
IF ( XINC  . N£ .  YINC) 

$  THEN 

PRINT  976 
STOP 
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M  M  U  M  W  M 


i 


i 


i 


t 


its. 

114. 

115. 

1 16  . 

117. 

1  IB  . 

119. 

120.  12 
121  . 

122  . 

123.  13 

124  . 

125. 

126. 

127  . 

128. 

129  . 

130.  15 

131  .  C 
132 .  C 

133  . 

134  . 

135. 

136. 

137  . 

138  . 

139  . 

140  . 

141  .  18 

142.  C 

143 .  C 
144  . 

145. 

146  . 

147  . 

148  . 

149  . 

150  . 

151  .  C 

152  .  C 

153  .  C 

154 .  C 

155 .  C 

156  . 

157  . 

158.  C 

159. 

160  . 

161  . 

162  . 

163  . 

164  . 

165  . 

166. 

167  . 

168. 

169. 


END  IF 

YOVbRX  «  SI2EY/S12EX 

XI  x  X0-O.5E0* (SI2EX-XINC) 

Y 1  =  Y0+0.5£0*(SI2£Y-Y1NC) 

RADIUS  *  SQRI( ( XINC*Y INC)/ PI ) 

XX( 1 )  *  XI 

DO  12  I =2, NUMX 

XXII)  =  XXI  1-1  )  +  XINC 

YY  (  1  )  =  Y  1 

DO  13  J  =  2 , NUMY 

YY  (  J  |  =  Y Y ( U- 1  ) "YINC 

M  =  0 

DO  15  J=1 , NUMY 
00  15  1=1. NUMX 
M  =  M+1 
X ( M  J  =  XX( I ) 

Y(f,1)  =  YY(J) 

CON  r 1MUE 

PRINT  SCHEMATIC  OF  DISTURBED  REGION  GRID 
I F ( NUMY  .GT.  10)  PRINT  915 
PRINT  920 

PRINT  930, (LABE13,J*1 , NUMX) 

00  18  K= 1 , NUMY 

PRINT  91 0 , LABELO  .  (  LABE  LI ,  J=1 .NUMX) 

PRINT  94  0, LABE  LO, ( ( K- 1 ) »NUMX+U , LABELO , J  =  1 ,NUMX) 

PRINT  910, LABELO, (LABEL1 ,  J  =  1 .NUMX) 

PRINT  9l0, LABELO, (LABEL2,J=1 ,NUMX) 

CONTINUE 

PRINT  COORDINATES  OF  GRID 
PRINT  920 
N  =  1 

WR I T  £ ( 6 , 100)  N,X(N),Y(N) 

WRITE (6, 100)  NUMX, X( NUMX) ,Y(NUMX) 

N  =  NU-NUMX+1 

WRITEI6, 100)  N , X ( N) , Y ( N) 

WR I T E ( 6 , 100)  NU.X(NU) , Y(NU) 

READ  GRID  DATA 

MODE  PARAMETERS  ARE  FROM  USING  NPUNCH=1  IN 
THE  WAVEGUID  OR  MODESRCH  COMPUTER  PROGRAMS 
READ  1000.10 
PRINT  1001,10 

AMD  I  ENT 

READ  1010, FREQ, SIGMA, EPSR 
PRINT  1 Q 1 1 .FREQ, SIGMA, EPSR 
OMEGA =6. 2831 85 3071 7959E3* FREQ 
WAVENO  =  .0209SB445EO*FREQ 
NGSQ  =  SIGMA/ ( IM»OMEGA*EPSO)+EPSR 
NGXTM  =  C  SQRT (NGSQ) 

IF(OREAL(NGXTM)  . LT .  O.OEO)  NGXTM*-NGXTM 
EXO  =  1 . OEO/NGXTM 

ECGNST *  ,03248E6*DM*EX0*WAVENO**2/(5.OE3*  SQRT (FREQ) ) 
N  =  0 


i 


i 
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170. 

PRINT  1009 

171  . 

READ  1020.THETA0.TMP1 

172. 

THTAO  =  THETA0*OTR 

173  . 

SO  =  C  SlN(THTAO) 

174. 

CSQ  =  (1 . EO , 0 . EO  )  -SO*  *2 

1 75  . 

C  =  C  SQRT(CSQ) 

176. 

SQROOT  =  C  SQR  T ( NGSQ~CSQ) 

177  . 

XTRAO  *  -4.0*TMP1*SO 

1 78  . 

PRINT  1021 .N.THETAO, XTRAO 

179  . 

KA  *  WAVENO*  RA  0 1  US 

100. 

KSO  =  WAVENO 'SO 

181  . 

ARGO  =  KSO*  R AD  I  US 

182  . 

CALL  CBESJY( ARGO, 1 ,BJ,BY,0,0) 

183. 

HI  2  =  BJ-  1 M*BY 

184  . 

<J1  =  BJ 

185  . 

CALL  CBESJY(ARG0,2,BJ,BY,0,0) 

186. 

H22  =  8J-IM*BY 

187  . 

J2  *  BJ 

188. 

C 

189. 

C 

DISTURBED 

190  . 

25 

READ  1010, F 

191  . 

I F ( F  .NE.  0.) 

192. 

$  THEN 

193. 

N  «  N*1 

194  . 

I F ( N  .EO.  MAXNX) 

195. 

S  THEN 

196. 

PRINT  974 

197. 

STOP 

198. 

ENDIF 

199. 

REAO  1020, THE TAP, TMP1 

200  . 

THTAP  *  THETAP*OTR 

201  . 

SP  *  C  SINlTHTAP) 

202.* 

XTRAP*-4.0*TMP1*SP 

203. 

PRINT  1021 .N.THETAP.XTRAP 

204  . 

S I ( N )  x  SP 

205. 

XTRI(N)  *  XTRAP 

206. 

GO  TO  25 

207  . 

ELSE 

208. 

NMAX  •  N 

209. 

END  IF 

210  . 

C 

211  . 

IF( NMAX  . EO.  1 ) 

212  . 

S  THEN 

213. 

C 

UNIFORM  DISTURBANCE 

214. 

NMAX  x  2 

215. 

S I ( 2 )  =  SP 

216. 

XTRK2)  *  XTRAP 

217  . 

ENDIF 

210  . 

C 

219. 

C 

SET  UP  INTERPOLATION  ARRAYS 

220. 

XDINC  x  0 . 5E0 • S I 2EX/( NMAX-1 ) 

221  . 

DO  33  N* 1 , NMAX 

222  . 

XI ( N)  x  XDINC* ( N-1 ) 

223. 

33 

Y I ( N )  *  YOVERX  * X I ( N ) 

224  . 

C 

225. 

c 

FILL  GRID  OF  SO  AND  XTRD 

226. 

IF ( I GRID  .EQ.  0) 
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$  then 


1 

1 


227  . 

229  .  C 

229  .  C 

230  .  C 

231  .  C 


232  . 

233  . 

234  . 

235  . 

236  . 

237  .  $ 

238  . 

239  . 

240  . 

241  . 

242  . 


243 .  C 

244.  35 

245  . 

246  . 

247.  37 

248  . 

249. 

250.  38 

251  .  S 


252  . 

253  . 

254  . 

255  . 

256  . 

257  . 

258  . 
259. 
260  . 


261  .  C 

262 .  39 

263 .  C 
264  . 

265 .  C 
266  . 

267  . 

268  . 

269 .  C 

270  . 

271  . 

272 .  C 

273  . 

274  . 

275  . 

276  . 

277  . 


278  . 

279  . 

280  . 
281  . 
282. 
233. 


SQUARE  OR  RECTANGLE 

A  RECTANGLE  CAN  ONLY  BE  UNIFORMLY  DISTURBED 
A  SQUARE  CAN  HAVE  A  NON-UNIFORM  DISTURBANCE 
LABEL  =  1  SQUARE  OR  RECTANGULAR  DISTURBANCE  * 
00  39  M=1,NU 
XM  «  A8S(X(M)-X0) 

YM  •  ABS(Y(M)-YO) 

IF(XM  .GT.  . 5E0»SI ZEX  .OR.  YM  .GT.  .5E0+SIZEY) 
then 

SD(M)  =  so 
XTRD(M)  *  XT  RAO 

ELSE 

I  •=  1 
K  =  1 


I  F ( XM  .LE.  X  I ( 1  +  1 ) )  GO  TO  37 
I  =  1  +  1 
GO  TO  35 

I  F ( YM  .LE.  Y I ( K+1 ) )  GO  TO  38 
K  a  K+1 
GO  TO  37 
I F ( I  .GT.  K) 

THEN 

SLOPE  a  (XM-XI ( I )  )/(XI(  1  +  1  ) -X I ( I ) ) 

SD(M)  a  SI(I )  +  (SI( 1  +  1)  —  SI(I) ) *  SLOPE 
XTRD(M)  a  XTRI( I )+(XTRI ( 1+1 )-XTRI( I ) )»SLOPE 

ELSE 

SLOPE  a  (YM-YI(K) )/( YI( K+1 )-YI (K) ) 

SD(M)  «  SI(K)+(SI (K+1 )— S I( K ) ) ‘SLOPE 
XTRD(M)  •  XTRI (K)+(XTRI (K+1 ) -XTR I ( K ) ) ‘SLOPE 

ENDIF 

ENDIF 

CONTINUE 


ELSE 

ELLIPSE  OR  CIRCLE 

LABEL  a  'CIRCULAR  OR  ELLIPTICAL  DISTURBANCE  ' 

A I  MAX  =  XI(NMAX) 

DO  59  M=1, NU 

CUCRi'S  OF  CENTER  OF  GRID  RELATIVE  TO  CENTER  OF  DISTURBANCE 

XM  =  X(M)-X0 
YM  a  Y(M)-YO 

COORDS  OF  CORNERS  OF  GRID 
XG(  1  )  a  XM-HXINC 
X(i(  2  )  a  XM+HXINC 
XG ( 3 )  a  XG( 1 ) 

XG(4)  =  XG( 2 ) 

Y  G (  1  )  =  YM+HY INC 

YG ( 2  I  =  YG ( 1 ) 

YG(3)  a  YM-H^INC 
YG ( 4 )  a  YC( 3 ) 

RGMAX  a  -1 . 0E6 
RGMIN  a  1.0E6 
DU  45  I » 1  , 4 
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r 


2 

284  . 

C 

RANGE 

OF 

EACH  CORNER  FROM  CENTER  OF  DISTURBANCE 

3 

285. 

RG 

* 

SORT ( XG( I )»*2+(YG(I)/Y0VERX)»*2) 

3 

286. 

RGMAX 

B 

AMAX1 (RGMAX.RG) 

3 

287  . 

45 

RGMIN 

B 

AMIN1 (RGMIN, RG) 

2 

288  . 

I F( RGMIN 

.GE.  AIMAX) 

2 

289. 

$ 

THEN 

2 

290. 

C 

GRID  IS  COMPLETELY  OUTSIDE  OF  DISTURBANCE 

3 

291  . 

SD(M)  *  SO 

3 

292. 

XTRD(M)  *  XTRAO 

3 

293. 

ELSE 

3 

294  . 

I F( RGMAX  .LE.  AIMAX) 

3 

295. 

$ 

THEN 

3 

296. 

C 

GRID  IS  COMPLETELY  INSIDE  OF  DISTURBANCE 

4 

297. 

1  =  1 

4 

298. 

A I  =  SORT ( XM*  *  2+ ( YM/YOVERX ) *  *2 ) 

4 

299. 

50 

I F( A I  .LE-  *1(1  +  1  ))  GO  TO  51 

4 

300. 

I  =  1  +  1 

4 

301  . 

GO  TO  50 

4 

302. 

51 

SLOPE  a  ( A I -X I ( I  )  )/(XI ( 1  +  1 )-XI(  I )) 

4 

303. 

SD ( M )  a  si ( I )  +  ( SI  ( 1  +  1 ) -SI ('  ) )»SLOPE 

4 

304. 

XTRD(M)  a  XTRI( I )+(XTRI( 1+1 , -XTRI( I ))*SLOPE 

4 

305. 

ELSE 

4 

306. 

C 

grid  IS  ON  BORDER  OF  THE  DISTURBANCE 

4 

307  . 

C 

DIVIDE  GRID  INTO  16  SUBSOUARES 

4 

308. 

DELTAX  =  XINC/8. 

4 

309. 

DELTAY  a  Y INC/8 . 

4 

310. 

XG( 1 )  »  ( XM-DE LTAX»3. )  »  »2 

4 

311  . 

XG( 2 )  *  (XM-DELTAX  )  •  *2 

4 

312. 

XG{3)  «  (XM+DELTAX  )  *  »2 

4 

313. 

XG( 4 )  «  (XM+DELTAX. 3. )»*2 

4 

314. 

YG( 1 )  a  ( ( YM+ DELTAY* 3. )/Y0VERX)*»2 

4 

315. 

YG( 2 )  »  ((YM+DELTAY  )/Y0VERX)»*2 

4 

316. 

YG(  3 )  •  ((YM-DELTAY  )/Y0VERX)**2 

4 

317. 

YG(4)  a  ( ( YM-DELTAY*3. )/Y0VERX)**2 

4 

318. 

C 

COUNT  NUMBER  OF  SUBSQUS  WITHIN  DISTURBANCE 

4 

319. 

N  a  0 

4 

320. 

DO  55  1*1,4 

5 

321  . 

XGSQ  *  XG(  I  ) 

5 

322. 

DO  55  J  =  1  .4 

6 

323  . 

RG  =  SORT(XGSO+YG(J)) 

6 

324  . 

55 

I F( RG  .LE-  AIMAX)  N*N+1 

4 

325. 

SLOPE  »  N/16. 

4 

326  . 

SD(M)  »  SO+(SI (NMAX)-SO)*SLOPE 

4 

327  . 

XTRD(M)  a  XTRAO+(XTRI( NMAX )~XTR AO)* SLOPE 

4 

328. 

ENDIF 

3 

329. 

ENOIF 

2 

330  . 

39 

CONTINUE 

1 

331  . 

END  IF 

1 

332. 

C 

1 

333. 

C 

SET 

UP 

1  A 

334. 

DO 

55 

Mai ,NU 

1 

335. 

DO 

65 

N* 1 ,NU 

2 

336. 

65 

A(M.N) 

«  CAPA(SD(N),SQRT((X(M)-X(N))**2+(Y(M)-Y(N) )**2) ) 

2 

337. 

C 

338. 

70 

MGR  ID 

.  0 

339. 

DO 

79 

M>  1 ,NU 

1 

340. 

NGRID 

«  0 

« 
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341  . 

XM  =  X  (M ) 

342  . 

YM  =  Y(M) 

HX INC  .GE.  O.OEO  .OR 

343. 

I F  (  XM+  HXINC  .LT.  O.OEO  .OR.  XM- 

344  . 

$  YM+  HYINC  .LT.  O.OEO  .OR.  YM~ 

HYINC  .GE.  O.OEO)  GO 

345  . 

C 

XM 1 R  13  INSIDE  GRID 

346  . 

C 

GRID  POINT  IS  DISTURBED 

347  . 

XG( 1  )  =  XM- . 25  *  X  INC 

348  . 

XG( 3 )  *  XG( 1 ) 

349  . 

XGI2)  =  XM* . 25  *  X I NC 

350  . 

XG ( 4 )  =  XG(2) 

351  . 

YGl  1  )  =  YM  * . 25  *  Y INC 

352. 

YG ( 2 )  =  YG( 1 ) 

353  . 

YG(3)  =  YM- . 25 *  Y INC 

354  . 

YG ( 4 )  r  YG(3) 

355. 

MGR  ID  =  M 

356  . 

PRINT  955, M 

357. 

NQUAD  =  0 

358  . 

SUM  =  0.010 

359. 

71 

NQUAD  =  NQUAD* 1 

360  . 

I F ( XG(NQUAD)* . 1 25*X1NC 

■  GE. 

O.OEO 

•  AND. 

361  . 

$  XGINQUADI-. 1 25*X INC 

.  LT  . 

O.OEO 

.AND  . 

362  . 

$  YG(NQUAD)*.  1 25  *  Y INC 

•  GE. 

O.OEO 

.AND. 

363  . 

S  YG(NQUAD|-. 125*YINC 

•  LT. 

O.OEO) 

GO  TO  71 

364  . 

NGR10  =  NCR  I D* 1 

365  . 

XM  =  XG( NQUAD) 

366  . 

YM  =  YGl NQUAD) 

367  . 

75 

RHOMN  =  SORT ( XM**2+YM**2) 

368  . 

ARGU  =  KSO  *  RHOMN 

369  . 

CALL  CBESJY(ARGU,1 .8U.BY.0, 

0) 

370  . 

HI  2  =  BJ” 1 M • 6Y 

371  . 

I F  (  XM  .EQ.  O.OEO)  XM-1  . 

0E-6 

372  . 

EZO(M)  =  C  SQR  T  <  ARGU*CAPG2( RHOMN ) ) * CON ST  »H12*XM/RHQMN 

373. 

I F ( NGRI D  .EQ.  0)  GO  TO 

78 

374  . 

SUM  -  SUM+EiO ( M ) 

375. 

IF( NQUAD  . LT .  4 )  GO  TO 

71 

376. 

EZO(M)  *  SUM/NGRID 

377  . 

78 

IF(C  ABS(EZOIM))  .LE.  1 

•  OE- 

6)  EZO ( M) = 1 .OE-21 

378  . 

79 

CONTINUE 

379. 

C 

380  . 

C 

PRINT  TABLE  OF  EZO 

3P1  . 

WRITE  (6, 1 10)  YMID 

382  . 

WRITE  (6.113) 

383  . 

DO  81  1=1, NUMY 

384  . 

WR I T  E ( 6 , 1  1 2 )  (C  ABS(EZO(NUMX*( 1-1 )*J)) ,J  =  1 .NUMX) 

385  . 

81 

CONT INUE 

386  . 

C 

387  . 

C 

SOLVE  TOR  SCATTERED  FIELDS 

388  . 

CALL  CL1NEQI A. EZO.EZS, IROW, 

CLNQ.NU.NDIM.NPTS, ERR) 

309  . 

C 

390  . 

( 

PRINT  TABLE  OF  EZS 

391  . 

WRITE  (6,111) 

392  . 

DO  83  1  =  1  .N'JMY 

393  . 

WR I T  E ( 6 . 1 12)  (C  ABS(EZS(NUMX*(I-1  )+J) ) . J*1 ,NUMX) 

394  . 

83 

CONT INUE 

395  . 

C 

396. 

DO  85  M= 1 ,NU 

397  . 

RATIO  .  EZS(M) /EZO(M) 
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398. 

399. 

400  . 

401  . 

402  . 
403. 
404  . 

405. 

406. 
407  . 

408. 

409. 

410  . 

411  . 

412. 

413. 

414. 
415  . 
416. 

417  . 

418  . 

419  . 
420. 

421  . 

422  . 

423  . 

424  . 

425. 

426. 

427  . 

428  . 

429. 

430. 
431  . 

432. 

433. 

434. 

435. 

436. 

437  . 

438  . 
439. 

440  . 

441  . 

442  . 
443. 

444  . 

445  . 
446. 
447  . 

448. 

449. 

450. 
451  . 

452. 

453. 

454. 


I F  (  MGR  I D  .EQ.  Ml  RATIO*1.0 
WI(M)  =  C  LOG(RATIO) 

CONTINUE 

PRINT  TABLE  OF  WI 
WR1 TE ( 6, 102) 

DO  87  1=1, NUMY 

MR  I T  E (6 , 103)  (REAL(MI<NUMX«(I-1)+U))»8.686,J*1 .NUMX) 

CONTINUE 

MRITE  (6,950) 

DIST  =  DMIN 
NPTS  =  NPTS+1 
I F ( NP  TS  .GT.  NRPTS) 

$  THEN 

PRINT  973 
STOP 

END  IF 

ARGU  *  KS0*01ST 
CALL  CBESJY(ARGU,1 .BJ.BY.0,0) 

HI  2  =  BJ-IWBY 

E2U  =  C  SQR  T ( A RGU*CAPG2( DIST)) +C0N5T  *H 1 2 

IF(XX(1 1-HXINC  .GE.  DIST  .OR.  XX ( NUMX ) +HX1 NC  .LE.  DIST  .OR. 

$  YY(1)+HYINC  .LE.  O.EO  .OR.  YY ( NUMY ) -HY I NC  .GE-  O.EO) 

$  THEN 

RCVR  IS  OUTSIDE  OF  DISTURBED  AREA 
XTRA2  =  XTRAO 
SUM  =  (O.OEO.O.OEO) 

DO  95  N=1 ,NU 

RHOMN  =  SORT (  (OIST-X(N) )**2+Y(N)**2) 

SUM  »  SUM+CAPA ( SD( N ) , RHOMN ) *E2S ( N) 

CONTINUE 
E2  *  EZU-SUM 
RATIO  =  EZ/EZU 

WMAG  *  20 . OEO * AL0G1 0 ( C  ABS(RATIO)) 

WANG  *  C  ANG( RATIO) 

ELSE 

RCVR  IS  IN  DISTURBED  AREA 
I  =  1 

DO  96  U=1 , NUMX  _ _ . 

IF(XX( J)-HXINC  .LT.  DIST  .AND.  DIST  .LE.  XX(d)+HXINC> 

$  GO  TO  97 

I  *  1  +  1 
CONTINUE 
K  *  I 

DO  98  U  =  1 , NUMY 

IF(YY( JJ-HYINC  .LT.  O.EO  -AND.  O.EO  .LE.  YY(U)+HYIN 
$  GO  TO  99 

K  *  K+  1 
CONTINUE 
M  ■  (K-1 )*NUMX+I 
XTRA2  =  XTRO(M) 

E2  »  EZU 

I F (M  , EQ-  MGRID)  GO  TO  130 
1 F( I  , GT .  1 )  GO  TO  120 
N1  ■  M 


.LT.  O.EO 


O.EO  .LE.  Y Y ( d ) +HY I NC ) 


MMIOMMUUUMMMMIOMUUUUUUUUWM 


105 


1 


455. 

456. 
457  . 

458. 

459. 

460  . 

461  . 

462  . 
463. 
464  . 

465. 

466. 

467  . 

468  . 
469. 

470  . 

471  . 

472  . 
473. 

474  . 

475  . 

476  . 

477  . 

478  . 
479. 

480  . 

481  . 

482  . 
483. 
484  . 

485. 

486. 
487  . 

488. 

489. 

490. 

491  . 

492  . 
493. 

494  . 

495  . 

496. 

497. 
498  . 
499. 

500  . 

501  . 

502  . 

503. 

504. 

505. 

506. 

507. 

508. 

509. 

510. 
511  . 


106 

107 

108 

109 

120 

121 

122 

123 

125 


$ 

S 

$ 


126 


N2  •  M  +  1 

I  F  ( K  .GT.  1 )  GO  TO  106 
N3  »  N1+NUMX 
N4  »  N  3+ 1 
GO  TO  125 

1 F ( K  .EQ.  NUMY)  GO  TO  108 
I F ( YY ( K)  .GE.  0. )  GO  TO  109 
N3  *  N1-NUMX 
N4  *  N3+1 
GO  TO  125 
N3  <  N1+NUMX 
N4  =  N  3+ 1 
GO  TO  125 

I F ( I  . LT.  NUMX)  GO  TO  122 
Ml  «  M-1 
N2  *  M 

IF(K  . EO.  1 )  GO  TO  105 
I F ( K  . EO.  NUMY)  GO  TO  108 
GO  TO  107 

I F ( XX (  I  )  .LE.  DIST)  GO  TO  123 

N1  «  M-1 

N2  «  M 

GO  TO  121 

N1  x  M 

N2  x  M+1 

GO  TO  121 

SLOPE  =  -Y(N1 )/( Y(N3)-Y(N1 ) ) 

Ml  «  MI (N1 )  +  (WI <N3)-WI ( N 1 ) )»SL0PE 
M2  «  MI(N2)+(WI(N4)-MI(N2) )«SLOPE 
I F ( 1  .EO.  NUMX  .AND.  DIST  .GT.  X ( N2 )  .OR. 

I  .EO.  1  -AND.  DIST  .LT.  X(N1)) 

THEN 

I F ( I  -EQ.  1) 

THEN 

DST  x  X0~.5*SIZEX 
DELTAX  x  HXINC 
M2  =  Ml 
N2  x  N 1 

ELSE 

DST  x  X0+ . 5*SI ZEX 
DELTAX  x  -HXINC 

END  IF 

ARGU  =  KSO  *  DST 

CALL  CBESJY( ARGU,  1  ,BJ . BY ,0,0) 

HI  2  =  BJ-IM*  BY 

EZ1  =  C  SORT ( ARGU* CAPG2 ( DST ) ) ♦ CONST »H1 2 
SUM  =  (O.CEO.O.OEO) 

DO  126  Nxj ,NU 

RHOMN  x  SORT ( ( DST  —  X ( N )  )  **2+Y ( N ) • * 2 ) 

SUM  x  SUM+C AP A ( SD ( N ) , RHOMN ) »EZS ( N ) 
CONTINUE 

Ml  «  CLOG( (1 .0,0.0 )-SUM/EZ1 ) 

RATIO  =  M1+(M2-M1 )*( DI ST-DST ) /DELTAX 

ELSE 

RATIO  «  Ml  +  ( M2~M 1 ) • ( (D I  ST— X(N1 ) ) /X  INC) 

ENOIF 

MMAG  «  REAL(RATI0)*8. 686 
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1 

1 

1 

1 

t 


1 

t 

1 

1 

1 

1 


1 

1 


512. 
513  . 

514. 

515. 

516. 

517. 

518. 

519. 

520  . 

521  . 

522. 

523. 

524. 

525. 

526. 

527. 

528. 

529. 
530- 
531  . 
532. 

533  . 

534  . 

535. 

536. 

537  . 

538  . 

539. 

540. 
541  . 

542. 

543. 

544. 

545. 

546. 

547. 

548. 

549. 

550. 

551  . 

552  . 

553  . 

554. 

555. 

556  . 

557  . 

558  . 
559. 

560  . 

561  . 

562. 

563. 

564. 

565. 

566. 
567  . 
568. 


WANG  «  AIMAG< RATIO) 

£2  »  C  £XP( RAT I0)*EZU 

EN01F 

C 

C 


IF(MGRID 

,  EQ. 

0) 

THEN 

XMTfi 

IS 

OUTSIDE  OF  DISTURBED  AREA 

XTRA1 

* 

XTRAO 

ELSE 

XMTR 

IS 

IN  DISTURBED  AREA 

XTRA1 

s 

XTRD(MGRID) 

END  IF 

C 

EZUABS  =£CONST*EZU«XTRAO 

EZPABS  s=£C0N5T  *EZ«C  SORT  ( XTRA1  *XTRA2) 

RATIO  .  EZPA8S/EZUASS 

EZREL  =  20 .OEO  » AL0G10(C  ABS( RATIO) ) 

EZANG  =  C  ANG( RATIO) 

EZUMAG  «  20.0E0»AL0G10(C  ABS(EZUABS)) 

EZUANG  *  C  ANG( EZUABS) 

EZPMAG  =  20 . OEO* AL0G1 0 (C  ABS(EZPABS)) 

E2PANG  «  C  ANfl ( E2PASS ) 

WRITE  16,960)  YMIO , 01  ST ,WMAG , WANG , EZRE L , EZANG , 

$  EZUMAG, EZUANG. EZPMAG, EZPANG 

YIPLOT(NPTS)  *  WMAG 
Y2PL0T(NPTS|  *  EZUMAG 
Y3P10T(NPTS)  •  EZPMAG 
I F ( I F LAG  .EQ.  2)  GO  TO  140 
C 

C  Y  VARIATION 

00  131  N  «  1 ,NU 

131  Y(N)  =  Y ( N ) +DE  LY 
DO  132  J-1.NUMY 

132  YY( J I  *  Y  Y  ( U J  +  DE LY 
XPLOT(NPTS)  *  YMID 
XLABEL  =  '  Y (KM)  ' 

YMID  =  YMUj+DEEY 

I F ( YMID  .LE.  YMAX)  GO  TO  70 
XMAX  =  YMAX 

IF(  I  PLOT  • EQ .  1 )  GO  TO  800 
STOP 
C 

C  DISTANCE  VARIATION 

140  XPLOT(NPTS)  =  DIST 

XLABEL  =  1  X (KM)  1 
OIST  =  0 1  ST  +  OELD 
IF ( DIST  . LE •  DMAX)  GO  TO  90 
XMAX  *  DMAX 

I F ( l PLOT  . EQ.  1 )  GO  TO  BOO 
STOP 
C 

C  PLOTTING 

800  CALL  DBP LOT 

C 

100  FORMAT ( '  COORDINATES  AT  CENTER  OF  MESH  NUMBER  ',13,'  ARE 

$  F8.2, '  Y«',F8.a) 
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569.  102  FORMAT !/,10X,'WI') 

570.  103  FORMAT(1X,13F10.4) 

571.  110  FORMA T ( 1  HO , 'YMID  «  '.F12.2) 

572.  Ill  FORMAT (/, 1 0X. * EZS' ) 

573.  112  FORMAT! ( IX, ( 1P13E10.2) )) 

574.  113  FORMAT!/, 10X, ' EZO' ) 

575.  910  FORMAT (20X,At , 13 AS) 

576.  915  FORMAT ('1') 

577.  920  FORMAT! 'O’) 

578.  930  F0RMAT(21X, 13A5) 

579.  940  FORMAT  (20X,A1 ,13(13, IX, A1)) 

580.  950  FORMAT!/ ,4X, 'YMID' ,5X. 'DIST 1 ,6X, '  WMAG  ',8X,'  WANG  ',8X,'  EZREL', 

581.  $  8X ,  '  EZ ANG 1 , 8 X ,  1 EZUMAG ‘ , 8X ,  1 EZUANG 1 , 6X ,  1 EZPMAG 1 , 8X ,  1 EZPANG 1 ) 

582.  955  FORMAT! 1 OXMTR  IS  INSIDE  GRID  ',13) 

583.  960  F0RMAT(2< 2X.F7 . 1  )  ,8(2X, 1PE12.5) ) 

584.  971  FORMAT!'  NU  *  NUMX«NUMY  IS  GREATER  THAN  PARAMETER  VARIABLE  NDIM 1 ) 

585.  972  FORMAT!'  NMAX  DOES  NOT  EQUAL  NU  OR  1') 

586.  973  FORMAT!'  NUMUER  OF  POINTS  PLOTTED  IS  GREATER  THAN  PARAMETER  VAR  I AB 

587 .  SLE  NRPTS ' ) 

583.  974  FORMAT! '  N  IS  GREATER  THAN  PARAMETER  VARIABLE  MAXNX ' ) 

589.  976  FORMAT!'  X1NC  MUST  EQUAL  YINC') 

590.  1000  FORMAT ( 20 A4 ) 

591.  1001  FORMAT! ' 1 ' ,20A4) 

692.  1009  F0RMAT(2X,  ' GR I  O'  ,6X, ' THETA'  ,25X,  '  XTRA'  ) 

593.  1010  FORMAT(10X,EB. 0 , 34X , E 1 0 . 0 , 2X , E5 . 0 ) 

594.  1011  FORMAT!'  FREQ  =  ',£10.3,'  SIGMA  «  *,E10.3,'  EPSR«',F7.2) 

595.  1020  FORMAT (1X,219.0,1X,2E15.0//) 

596.  1021  FORMAT (1X,13,2F10.5,2XI1P2E18.9) 

597 .  END 

END  FTN  2678  IBANK  22972  DBANK  2068  COMMON 
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SUBROUTINE  DliPLOT 
PARAMETER  NRPTS=500 
CHARACTER*  8  XLA8CL 
CHARACT ER* 36  LABEL 
CHAR ACT  ER*  4  10(20) 

01  MENS  I  ON  XPLOT(NRPTS) ,71 PLOT(NRPTS) , Y2PL0T ( NRPTS ) , Y3PL0T ( NRPTS) 
COMMON/ LABELS/ID, LABEL, XLABEL 

COMMON/ INPUT/ FREQ, OMIN.OMAX.SIZ EX. SI ZEY.XO, YO , EXT  I C , E> T IC , WXT IC , 
$  WYT I C , XLNG, Y  LNG, EMIN, EM AX , WMIN , WMAX , XMAX , NUMX , NUMY , 

S  I  FLAG 

COMMON/ PLOATA/XPLOT , Y1PL0T ..Y2PL0T, Y3PL0T.NPTS 
CALL  BGNPLd  ) 

CALL  PHYS0R(2. 0,2.0) 

CALL  TITLE('  1 , 1 .XLABEL, B, 'WMAG(DB) 1 ,8, XLNG, YLNG) 

CALL  INTAXS 
CALL  YAXANG(O.O) 

CALL  GR A F ( 0 . .WXTIC, XMAX, WMIN, WYTIC, WMAX) 

CALL  CL'RVElXPLOT.YIPLOT.NPTS.O) 

CALL  MESSAG( 10,80,0.0,-0.8) 

CALL  MESSAG( LABEL, 36, 0. 0,-1 .0) 

CALL  MESSAG( 1  FREQ  =  ' , 7 , 0 . 0 , - 1 . 2 ) 

CALL  REALNO(  FREO, 3,0. 7,-1 ,2) 

IF(IFLAG  .EQ.  1)  CALL  MESSAGf 'OMIN  « 1 , 6 , 2 . 0 , -1 . 2 ) 

IF(IFlAG  .EQ.  1)  CALL  REALNO(  DMIN , 1 , 2 . 7 , -1 . 2 ) 

CALL  ME SSAG( ' S I ZEX -  X0« 

CALL  REALNO(  SIZEX, 1 ,0.7 ,-1 .4) 

CALL  REALNO(  XO , 1 , 2 . 0 , - 1 . 4 ) 

CALL  INTN0(NUMX,3.7,-1 .4) 

CALL  MESSAG( ' SIZEY*  Y0» 

CALL  RCALNO(  S I ZEY , 1 . 0 . 7 , - 1 . 6) 

CALL  REAlNO(  Y0,1 ,2. 0,-1 .6) 

CALL  INTNO(NUMY,3.7,-1 .6) 

CALL  ENDPL( 1 ) 

CALL  BGNPL ( 2 ) 

CALL  PHYS0R(2. 0,2.0) 

CALL  TITLE!'  ' , 1 .XLABEL. 8. 'OB  ABOVE  1  MI CROVOLT/ME TER ' , 26 , 

S  XLNG, YLNG) 

CALL  INTAXS 
CALL  YAXANG(O.O) 

CALL  GRAF(0. , EXT IC, XMAX, EMIN. EYTIC.EMAX) 

CALL  CURVE (XPL0T,Y2PL0T.NOTS,0) 
call  dash 

call  CURVE(XPLOT, Y3PLOT,NPTS,0) 
call  R2SET( 'DASH' ) 

CALL  STRTPT(XLNG-2.0,YLNG) 

CALL  CONNPT(XLNG-0.9,YLNG) 


call  MESSAGl 1  EZUMAG' , 16.XLNG-2.0, YLNG) 

CALL  MESSAG( ' -  EZPMAG’, 1 6 , X LNG-2 . 0 , Y LNG-0 . 3 ) 

CALL  MESSAG( 10,80 ,0 . 0,-0 . B) 
call  MESSAG( LABEL, 36,0 .0 .-1 .0) 
call  MESSAG( ' FREO  =  ' , 7 , 0 . 0 , - 1 . 2 ) 

CALL  REALNO(  FREO, 3.0. 7,-1 .2) 


1F(IFLAG  .EQ.  1)  CALL  M£SSAG('DMIN  » ' , 6 , 2 . 0 , - 1 . 2 ) 
JFIIFLAG  .EQ.  1)  CALL  RE A LN0(  OMIN , 1 , 2 . 7 , -1  . 2 ) 


NUMX*  ' ,32.0. 0,-1 .4) 

NUMY=  * ,32,0. 0,-1 ,6) 
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55 

CALL 

MESSAG( 1 SI2EX-  XC» 

‘1UMX- 

1 ,32 ,0. 0,-1 .4) 

56 

CALL 

REALNO(  SIZEX, 1 ,0.7 ,-1 .4) 

57 

CALL 

REALNO(  XO.1,2.0,-1 .4) 

58 

CALL 

INTNO(NUMX,3. 7,-1 .4) 

59 

CALL 

MESS AG ( ' SIZEY=  ¥0= 

NUM  Y = 

* ,32,0. 0,-1 .6) 

60 

CALL 

REALNO(  SIZEY.1 ,0.7, -1.6) 

61 

CALL 

REA  LNO(  YO, 1 ,2. 0,-1 .6) 

62 

CALL 

INTNO(NJMY,3.7,-1 .6) 

63 

CALL 

ENDP  L ( 2 ) 

64  RETURN 

65  END 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 


SUBROUTINE  CBESJY ( 2 , K , SJ , BY , KIND , NPR I  NT ) 

IMPLICIT  COMPLEX  (A-H.O-Z) 

REAL  DGMSUM.DG1 , DG2.P 1/3.1 4159265358979E0/ ,RRT , 

S  EULER/O. 577215664901533E0/ 

THIS  SUBROUTINE  CALCULATES  BESSEL  FUNCTIONS  OF  THE  FIRST  KIND(JN)  AND 
OF  THE  SECOND  KINO(YN).  N  IS  THE  ORDER  AND  IS  ALLOWED  TO  BE  ANY 
POSITIVE  INTEGER. 

THE  ARGUMENT  MUST  BE  DECLARED  COMPLEX  IN  THE  CALLING  ROUTINE. 

IF  A  REAL  ARGUMENT  IS  DESIRED  JUST  SET  THE  IMAGINARY  PART  TO  ZERO. 

FOR  THE  ARGUMENT  Z«RHO*C  EXP(I»PH1)  THE  EMPLOYED  EQUATIONS  ARE  VALID 
AS  FOLLOWS: 

IF  0<RH0< 1 3  THEN  PH  I  =  ANY  VALUE 

IF  RHO* 1 3  OR  1 3<RH0< INF  I NI T Y  THEN  PHI*PI  OR  -PI<PHI<PI 

FOR  RHO=0  JO  AND  J1  ARE  SET  TO  THEIR  CORRECT  VALUES. 

HOWEVER,  YO  AND  Y1  ARE  NOT  CALCULATED  SINCE  THEY  APPROACH  INFINITY 
IN  THE  REAL  AND/OR  IMAGINARY  PART.  THEREFORE  AN  ERROR  MESSAGE  IS 
PRINTED. 


Z  IS  THE  ARGUMENT. 

K  IS  THE  ORDER. 

BJ  IS  THE  BESSEL  FUNCTION  OF  THE  FIRST  KIND. 

BY  IS  THE  BESSEL  FUNCTION  OF  THE  SECOND  KIND. 

KIND* 1  CAUSES  CALCULATION  OF  THE  BESSEL  FUNCTION  OF  THE  FIRST  KIND 
ONLY  . 

ANY  OTHER  VALUE  OF  KIND  CAUSES  CALCULATIONS  TO  BE  DONE  FOR  BESSEL 
FUNCTIONS  OFF  BOTH  THE  FIRST  AND  SECOND  KIND. 

NPRINT=0  CAUSES  NO  DEBUG  PRINTOUT. 

NPR i NT* 1  CAUSES  DEBUG  PRINTOUT. 

THE  CALLING  STATEMENT  MUST  HAVE  DECLARED  THE  PARAMETERS  CORRESPONDING 
TO  Z,  BJ,  AND  BY  AS  COMPLEX 


I F ( C  ABS(Z)  . NE .  O.EO)  GO  TO  7 
B J* ( 0 . OEO , 0 . OEO) 

I F ( K  .EQ.  0)  BJ«(1 .OEO.O . OEO) 

I F ( KIND  .NE.  1  .AND.  K  .EQ.  0)  PRINT  400 
400  FORMAT! 1H0, •**«  Y  NOT  CALCULATED  FOR  ARGUMENT  OF  MAGNITUDE  0’//) 
RETURN 

7  I F ( C  ABS(Z)  .LT.  13. EO)  GO  TO  10 


ASSYMPTOTIC  EXPANSION 

RH0*8 , »  Z 
MU*4»K««2 

RT*C  SQRT ( 2. /( PI *Z ) ) 
RRT=RT 

I F ( RRT  .LT.  0.)  RT*-RT 
P«0. 

DO  LOOP  FOR  CALCULATING  P 
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55 

00  1  N= 1 , 30 

56 

M  =  N-1 

57 

MM=2*M 

56 

I F ( N  .EO.  1}  GO  TO  ? 

59 

TERM=(-1 )*(MU-(2*MM-3)*»2)«(MU-(2' 

60 

(TERM 

61 

I F ( NPRI NT  .EO.  1)  PRINT  100, TERM 

62 

100 

FORMAT ( IX, 'TERM-' .2E30.15) 

63 

P.P+TERM 

64 

I F ( C  ABS(TERM)  .GE.  C  ABS(TERMS) 

65 

$TO  3 

66 

TERMS*  TERM 

67 

GO  TO  1 

68 

2 

T  ERM= 1 . EO 

69 

P*TERM 

70 

1 F ( NPR I NT  .EO.  1)  PRINT  100, TERM 

71 

T  ERMS*T  ERM 

72 

1 

CONTINUE 

73 

3 

CONTINUE 

74 

I F ( NPR I  NT  .EO.  1)  PRINT  200, P 

75 

200 

FORMAT! 1H0,'P*‘ .2E30.15) 

76 

0  =  0. 

77 

C 

DO  1 

LOOP  FOR  CALCULATING  0 

78 

DO  4  N=1 ,30 

79 

M  =  N-1 

80 

MM=2*M 

81 

MMM=2*M+ 1 

82 

IF(N  .EO.  1 )  GO  TO  5 

83 

TERM* ( - 1 )*(MU-(2*MM-1 )**2)»(MU-(2 

84 

STERM 

85 

I F ( NPRI NT  .EO.  1)  PRINT  100, TERM 

86 

0=0+TERM 

87 

I F ( C  ABS (TERM )  .GE.  C  ABS( TERMS ) 

88 

$  TO  6 

89 

TERMSsTERM 

90 

GO  TO  4 

91 

5 

TERM* (MU-1 )/RHO 

92 

Q=TERM 

93 

1 F ( NPR I  NT  .EO.  1)  PRINT  100, TERM 

94 

TERMS-TERM 

95 

4 

CONTINUE 

96 

6 

CONTINUE 

97 

I F ( NPRI NT  .EO.  1)  PRINT  300,0 

98 

300 

FORMAT! 1H0, '0=' .2E30.15) 

99 

Bd=RT«C  C0S(Z-K*PI/2.-PI/4. )*P-RT 

100 

I F ( K I NO  .EO.  1 )  GO  TO  8 

101 

6Y*RT*C  SIN(Z-K*PI/2.-PI/4. J*P+RT 

102 

8 

RETURN 

103 

C 

104 

C 

105 

C 

POWER  SERIES  EXPANSION 

106 

c 

107 

10 

NTERMS=35 

108 

KFAC*1 

109 

I F  ( K  .LE.  1  )  GO  TO  30 

110 

00  20  N*2,K 

111 

20 

kfac«kfac«n 
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112 

30 

TERM* ( Z/2 . )**K/KFAC 

113 

BJ*  TERM 

114 

I F ( K IND  .EQ.  1 )  GO  TO  91 

115 

DG1  *0. 

116 

DG2*0. 

117 

I F ( K  .EQ.  0)  GO  TO  BO 

118 

DO  60  N=1 ,K 

119 

60 

DG2*DG2+ 1 ./N 

120 

80 

T  SUM3=_  TERM*  DG2 

121 

SUMT3*TSUM3 

122 

91 

DO  40  M= 1 .NTERMS 

123 

TERM* TERM* (1/2. )»*2*(-1 ) / { ( K+M) •«) 

124 

BJ*BJ+TERM 

125 

I  F  ( KINO  .EQ.  1 )  GO  TO  92 

126 

DG1 *DG1 + 1 .EO/M 

127 

DG2*0G2+ 1 . EO/(M+K) 

128 

DGMSUM=DG1+0G2 

129 

TSUM3*-TERM*0GMSUM 

130 

SUMT3=SUMT  3+TSUM3 

131 

92 

I  F ( C  ABS ( TERM )  . LE .  1 . E- 1 7 )  GO  TO  50 

132 

40 

CONTINUE 

133 

50 

I F ( KINO  .EQ.  1 )  GG  TO  93 

134 

TERM3=SUMT3/PI 

135 

T ERM1  =  ( 2  .  /  P I )*BJ*(£ULER+C  L0G(Z/2. )) 

136 

SUMT2= ( 0 . ,0. ) 

137 

I F ( K  -EQ.  0)  GO  TO  120 

138 

KM1FAC=KFAC/K 

139 

TSUM2*KM1FAC*(Z/2. )**(-k) 

140 

S'JMT2=TSUM2 

141 

I F ( K  -EQ.  1)  GO  TO  120 

142 

KM1 *K— 1 

143 

DO  130  Ms  1 ,KM1 

144 

KMM=K-M 

145 

TSUM2*TSUM2/( * (Z/2 . >**2 

146 

130 

SUMT2=S'JMT2+TSUM2 

147 

120 

TERM2.-SUMT2/PI 

148 

BY.TERM1 +TERM2+TERM3 

149 

93 

RETURN 

150 

END 
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1  SUBROUTINE  CLIN  EQ  ( A , B, X , 1  ROW, Q .N , NDIM, I  FLAG, ERR ) 

2  C 

3  C  CUN  EO  USES  L-U  DECOMPOSITION  TO 

4  C  FIND  THE  TRIANGULAR  MATRICES  L,  U 

5  C  SUCH  THAT  L  «  U  «  A.  L  AND  U  ARE 

6  C  STORED  IN  A.  THIS  FORM  IS  USED  WITH 

7  C  BACK-SUBSTITUTION  TO  FIND  THU  SOLN 

6  C  X  OF  A*X-L*U*X*B. 

9  C  N  IS  THE  NUMBER  OF  EQUATIONS  AND 


10 

C 

N 

DIM  IS  THE  DIMENSION  OF  ALL  ARRAYS 

11 

C 

IN 

THE  PARAMETER  LIST. 

12 

c 

13 

c 

IF 

I FLAG  *  0,  L,  U.  AND  X  ARE 

14 

c 

COMPUTED. 

15 

c 

IF 

I FLAG  IS  NON-ZERO.  IT  IS  ASSUMED 

16 

c 

THAT  L  AND  U  HAVE  BEEN  COMPUTED  IN 

17 

c 

A 

PREVIOUS  CALL  AND  ARE  STILL  STORED 

18 

c 

IN 

A.  THUS  ONLY  X  IS  COMPUTED. 

19 

c 

ERR  IS  THE  ESTIMATED  RELATIVE 

20 

c 

ERROR  OF  THE  SOLUTION  VECTOR. 

21 

c 

22 

COMPLEX  A.  8,  X.  T 

23 

REAL  ERR 

24 

DIMENSION  A(NDIM.NDIM) , B ( NDIM ) , X ( NDIM ) 

25 

DIMENSION  IROW(NDIM) ,0< NDIM) 

26 

DATA  EPS  /I . 0E“15/ 

27 

c 

28 

c 

29 

IF  (N.GT.NDIM)  GO  TO  900 

30 

IF  ( I FLAG. NE . 0 )  GO  TO  600 

31 

DO  050  1=1 , N 

32 

0(1)  »  0.0 

33 

DO  040  J  *  1 ,N 

34 

00  =  C  ABS  ( A (  I , J ) ) 

35 

040 

IF  (O(I).LT.OO)  Od)  *  00 

36 

IF  (0(1). EO. 0.0)  GO  TO  901 

37 

050 

CONTINUE 

38 

ERR  =  EPS 

39 

PPIV  =  0.0 

40 

DO  100  I  «  1  ,  N 

41 

100 

I  ROW ( I )  =  I 

42 

c 

43 

DO  500  L  *  1  .  N 

44 

PIVOT  *  0.0 

45 

K  ■  L  —  1 

46 

DO  240  I  =  L.N 

47 

IF  (K.LT.1)  GO  TO  230 

48 

DO  220  J  •  1 , K 

49 

220 

A(I.L)  =  A(I.L)  ~  A ( J , L )  *  A( I , d ) 

50 

230 

F  =  C  ABS  ( A( I , L) )  /  0(1) 

51 

IF  (PIVOT. GT.F)  GO  TO  240 

52 

PIVOT  «  F 

53 

NPIVOT  «  I 

54 

240 

CONTINUE 

54 


55 

IF  (PIVOT. £0.0.0)  GO  TO  901 

56 

IF  (PPIV.LE. PIVOT)  GO  TO  250 

57 

ERR  *  ERR  *  PPIV  /  PIVOT 

58 

IF  (ERR.GE.i .0)  GO  TO  901 

59 

250 

PPIV  x  PIVOT 

60 

IF  (NPlVOT.Eg.L)  GO  TO  280 

61 

O(NPIVOT)  *  Q(l) 

63 

J  •  I ROW! L) 

63 

IROW(L)  =  IROW(NPIVOT) 

64 

IROW(NPIVOT)  *  J 

65 

DO  260  I  •  1 ,N 

66 

T  =  A  ( l ,  I  ) 

67 

A ( L , I )  =  A(NPIVOT.I) 

68 

A(NPIVOT.l)  «  T 

69 

260 

CONTINUE 

70 

280 

IF  (L.EQ.N)  GO  TO  500 

71 

T  «  (1 .OEO.O.OEO)  /  A(L,L) 

72 

K  *  L  +  1 

73 

M  «  L  -  1 

74 

DO  450  I  *  K , N 

75 

IF  (M.I.T.1  )  GO  TO  400 

76 

00  350  J  -  1  ,M 

77 

350 

A ( L , I )  =  A(L.I)  -  A ( L<  0)  *  A(U.I) 

78 

400 

A(L,I)  *  T  *  A(L.I) 

79 

450 

CONTINUE 

80 

500 

CONTINUE 

81 

IF  (ERR.GT.1 .OE-5)  PRINT  998,  ERR 

82 

C 

83 

C 

84 

600 

DO  620  I  «  2,N 

85 

620 

X(I)  *  ( 0 . OEO ,0.060) 

86 

J  *  IROW( 1 ) 

87 

X ( 1 )  x  B(d)  /  A(1,1) 

80 

DO  700  I  »  2,N 

89 

J  x  IRQW(I) 

90 

K  «  I  -  1 

91 

DO  650  L  ■  1 , K 

92 

650 

X ( I )  »  X(I)  +  A(I,L)  *  X (  L) 

93 

X ( I )  x  (8(d)  -  X C I ) )  /  A ( I ,  I ) 

94 

700 

CONTINUE 

95 

K  =  N  -  1 

96 

DO  800  I  x  1  ,K 

97 

d  =  N  -  I 

98 

M  =  d  +  1 

99 

DO  800  L  •  M,N 

100 

X(d)  «  X ( d )  -  X(L)  *  A(d , L) 

101 

800 

CONTINUE 

102 

RETURN 

103 

C 

104 

900 

PRINT  999 

105 

ERR  *1.0 

106 

RETURN 

107 

901 

PRINT  997 

108 

ERR  »  1.0 

109 

RETURN 

HO 

997 

FORMAT  ('1  ERROR  IN  CLIN  EQ,  MATRIX 

111 

998 

FORMAT  ('  CAUTION-', 
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SINGULAR') 


$  1  CLIN  EQ  HAS  DECOMPOSED  AN  ILL-CONDITIONED  MATRIX.',/, 

$  '  RESULTS  WILL  HAVE  RELATIVE  ERROR  >',E11.2) 

999  FORMAT  ( 1 1  ERROR  IN  CLIN  EQ,  MATRIX  SIZE  GREATER  THAN  NDIM' ) 
END 
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(OOB«4aiUIAUK) 


r 


i 


1  FUNCTION  C  ANG(ARG) 

IMPLICIT  REAL  (A-H.O-Z) 

COMPLEX  ARG,MINUSI/(O.OEO,-1 .OEO)/ 
ARGR«ARG 
ARGI»MINUSI*ARG 
C  ANG*  ATAN2( ARGI ,ARGR) 

I F( ARGI  . LT.  o.)  C  ANG«C  ANG+6. 2831853072E00 
RETURN 
ENU 
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1 

COMPLEX 

FUNCTION  CAPA(SP.RHO) 

2 

IMPLICIT  COMPLEX  (A-H.O-Z) 

3 

COMPLEX 

IM/(0.£0. 1  . EO ) / , J1  ,<>3 

4 

COMPLEX 

KSO 

5 

REAL  KA 

■  RHO 

6 

COMMON/ONE/HI  2 ,H22 , J1 , J2 , SO , ARGO , KA , KSO 

7 

DEFINE 

CAPGI(ARG)  *  SORT ( ARG/( 6371 .OEO*  SIN(ARG/6371 .OEO) ) ) 

8 

ARGP  « 

KA*SP 

9 

TERM1  > 

( 0 . OEO ,0.7853981 63397448E0 ) * ( ( SP/SO ) **2-1 .OEO) 

10 

TERM2  - 

2. OEO* ARGO* ( 1 .OEO- . 25E0* ARGP*  *2 ) 

11 

I F ( RHO 

.EO.  0.0)  THEN 

12 

CAPA  •  ( SP/SO )**4*-TERM1*i TERM2*H12+ARGP**2*H22) 

13 

ELSE 

14 

COFAMN  =  TERM1  *  (  TERM2*  J1*-ARGP**2*J2 ) 

15 

ARGU  •  KSO* RHO 

16 

CALL  CBESJY(ARGU.O.BJ.BY.O) 

17 

H02  *  BJ-1M«BY 

18 

CAPA  •  COFAMN*H02»CAPG1 (RHO) 

19 

END  IF 

20 

RETURN 

21 

END 
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